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ABSTRACT 


A variation  of  Gaussian  elimination  Is  presented  for  solving  band 
systems  of  linear  equations  on  computers  with  limited  core  storage, 
without  the  use  of  auxiliary  storage  such  as  disk  or  tape.  The  method 
Is  based  on  the  somewhat  unusual  Idea  of  recomputing  rather  than  saving 
most  nonzero  coefficients  In  the  reduced  triangular  system,  thus  trading 
an  Increase  In  work  for  a decrease  In  storage.  For  a five-point  problem 


on  an  n X n grid,  the  storage  required  !•  “ n^  versus  “ for  band 
elimination,  while  surprisingly  the  work  required  at  most\  doubles. , 


k 

j 


f 
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1.  Introductton 


\ 


Consider  the  system  of  llnesr  equations 
(S)  A X - b 


where  the  coefficient  matrix  A Is  an  N x N symmetric  positive  definite 
band  matrix  with  bandwidth  m,  l.e.,  a^j  ■ 0 if  |1  - J|  > m (see  Figure 

1).  Direct  methods  for  solving  (S)  are  generally  variations  of 

symmetric  Gaussian  elimination:  Form  the  U^DU  decomposition  of  A,  where 
U Is  unit  upper  triangular  and  D positive  diagonal,  and  then 
successively  solve  the  triangular  systems 

U*’y“b,  Dz“y,  Ux“z. 


The  matrix  U is  again  a band  matrix  with  the  same  band%d.dth  as  A. 
Variations  of  symmetric  Gaussian  elimination  which  take  advantage  of 
this  structure  to  avoid  storing  and  operating  on  entries  outside  the 
band  are  known  as  band  elimination  methods  [6],  The  total  work  (in 
terms  of  the  number  of  multiplications  and  divisions)  and  storage 
required  are  given  by 

e_(N,m)  - 1/2  Nm^  + 7/2  Nm  - 1/3  m^  + 0(N+m^) 

O 

and 

S„(N,m)  - Nm  + 0(N+in^) 

B 

respectively  (IJ. 

As  an  example,  consider  the  following  model  problem  which  arises 
from  the  familiar  five-point  finite  difference  discretization  of  the 
Poisson  equation  on  the  unit  square  with  homogeneous  boundary 
conditions.  Given  a uniform  n x n grid  in  the  plane  (see  Figure  2a),  we 
associate  a variable  u^j  with  each  mesh-point  (i,J)  and  form  the  system 

of  linear  equations 


4 u 


ij 


■ Vi,j  ■ Vl,j 


■ “i.J-1  ■ “ij+l  " 


1^1,  j <,  n 


where 


Uij  =0 


i » 0,  n+1  or  J ” 0,  n+1 


n unknown  variables  , and,  with  the  natural  row-by-row 


There  are  N 

ordering  (see  Figure  2b),  the  bandwidth  is  m » n 

•storage  required  to  solve  the  linear  system  are 
soectlvely. 


Thus  the  work  and 
1 /2  n^  and  n^ 
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Our  model  problem  Illustrates  the  behavior  one  often  encounters  In 
using  Gaussian  elimination  to  solve  large  band  systems:  The  storage 
required  can  easily  exceed  the  core  storage  available  for  even 
moderately  large  N and  m,  even  though  the  problem  and  solution  (l.e..  A, 
b,  and  x)  CAN  be  represented  In  core.  Thus,  although  we  could  store  the 
nonzero  coefficients,  right  hand  side,  and  solution  for  our  model 
2 

problem  in  ' Sn  locations,  the  factorization  would  require  an 


additional  ' n locations.  More  generally,  although  we  could  store  x In 
N words  of  storage,  we  would  still  require  an  additional  ~ Nm  words  of 
storage  for  D and  the  band  of  U,  even  If  we  could  recompute  the  nonzero 
entries  of  A and  b as  easily  as  store  them. 


In  this  paper  we  discuss  several  variations  of  band  elimination 
which  can  solve  (S)  with  minimal  core  storage.  The  methods  are  based  on 
the  following  assumptions: 

(1)  The  nonzero  entries  of  A and  b are  Inexpensive  to  generate  on 
demand  (e.g.,  A is  sparse  and  can  be  represented  more  compactly 
than  in  band  form,  or  A must  be  preserved  for  subsequent 
computations! . 


(2)  There  Is  enough  core  storage  for  the  solution  vector  x,  plus  an 

2 

additional  yCN-Hn  ) words  of  working  storage  (for  some  small 
constant  y) , but  not  enough  to  store  the  band  of  U. 

We  count  only  the  working  storage  required  to  solve  (S)  In  stating  the 
storage  requirements  of  such  methods.  All  other  storage  (l.e.,  storage 
for  A,  b,  and  x)  Is  associated  with  the  linear  system  rather  than  the 
method  of  solution  and  is  Ignored. 


In  section  2,  we  review  the  standard  approach  to  this  problem,  the 
use  of  auxiliary  storage,  and  Indicate  some  of  the  disadvantages  of 
auxiliary  storage  band  elimination.  In  section  3,  we  Introduce  minimal 
storage  methods  as  a corollary  of  the  somewhat  unusual  Idea  of 
recomputing  rather  than  saving  most  nonzero  entries  In  the 

factorization,^  the  effect  being  to  trade  a logarithmic  Increase  In  work 
for  a factor  of  N/2m  decrease  in  storage.  Finally,  in  Section  4,  we 
consider  the  special  case  of  the  model  problem,  and  show  that  the  work 

3 2 

at  most  doubles  while  storage  Is  reduced  from  " n to  ~ n . 


The  same  approach  can  be  applied  to  sparse  elimination;  see  [2], 


J 
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2.  Auxiliary  Storage  Band  Elimination 

A standard  approach  to  problems  in  which  the  storage  required 
exceeds  the  amount  of  core  storage  available  Is  to  use  auxiliary  storage 
such  as  disk  or  tape.  This  approach  works  extremely  well  for  band 
elimination  (5],  as  can  be  seen  by  examining  the  process  more  closely. 

In  the  classic  view  of  Gaussian  elimination,  we  use  the  k^^ 

equation  to  eliminate  the  k^^  variable  from  the  remaining  N-k  equations 
for  k - 1,...,N,  and  then  solve  the  resulting  triangular  system.  In 
terms  of  the  factorization  model,  this  corresponds  to  solving 

U*’  y - b,  D - y 

as  A Is  factored,  and  then  solving 

U X - z 

when  the  factorization  Is  complete. 

A program  fragment  for  symmetric  band  elimination  using  a scratch 
array  M appears  as  Figure  3a.  Note  that  only  those  locations  M(1,J]  for 

which  k £ 1 j<  j £ k+m  are  referenced  in  eliminating  the  k*"^  variable. 
These  locations  form  an  (nrt-l  )x(iitf  1 ) triangular  window  on  the  upper  band 
of  M.  All  the  previous  rows  contain  entries  of  U and  will  not  be  needed 
until  back-solution;  all  subsequent  columns  contain  entries  of  A and 
have  not  yet  been  used  In  the  elimination  process;  and  the  window 
advances  as  each  variable  is  eliminated  (see  Figure  3b). 

This  suggests  an  auxiliary  storage  band  elimination  algorithm.  We 
use  core  storage  as  an  (mfl )x(nH-l ) triangular  window  on  the  elimination 
process.  After  each  variable  Is  eliminated,  the  window  Is  shifted;  the 
row  of  U left  behind  Is  written  to  auxiliary  storage;  and  the  last 
column  of  the  window  Is  Initialized  to  the  next  column  of  A.  During 
back-solution,  we  retrieve  the  rows  of  U in  reverse  order  and  solve  for 
the  corresponding  component  of  the  solution.  Since  we  do  exactly  the 
same  operations  as  In  band  elimination, 

® ^ Nm^, 

Ab 

yet  the  core  storage  required  (other  than  that  for  A,  b,  and  x which  we 
have  agreed  to  Ignore)  Is  just  that  for  the  window,  namely 

S^^(N,m)  - 1/2  (nH-l)(nri-2). 

AS 


An  obvious  refinement  of  this  algorithm  is  to  save  as  much  of  the 
factorization  In  core  as  possible,  but,  when  core  storage  Is  exhausted, 
to  make  room  by  moving  those  rows  of  U not  needed  for  the  next  stage  of 
elimination  (i.e.,  not  In  the  window)  to  auxiliary  storage.  In  a 
certain  sense.  In-core  band  elimination  on  computers  with  large  virtual 


memory  has  the  same  effect:  When  no  real  memory  la  available,  the 
operating  system  moves  the  least  recently  used  pages  (l.e.,  rows  of  U) 
out  to  the  swapping  device  (l.e.,  auxiliary  storage),  and  then  retrieves 
them  as  they  are  referenced  during  back-solution. 

While  auxiliary  storage  band  elimination  can  be  used  to  solve  large 
band  systems  of  linear  equations  on  computers  with  limited  core  storage. 
It  does  have  certain  disadvantages.  First,  a large  amount  of  auxiliary 
storage  Is  required  to  save  the  factorization.  Second,  the  added  cost 
of  transferring  the  factorization  to  and  from  auxiliary  storage  can  not 
be  neglected.  While  In  principle  this  Input/output  can  partially 
overlap  the  numerical  computation.  In  practice  the  degree  of  such 
overlap  is  highly  dependent  on  the  machine,  auxiliary  storage  device, 
and  operating  system  In  use.  Third,  retrieving  the  rows  of  U In  reverse 
order  requires  the  ability  to  read  backward  or  random  access  auxiliary 
storage,  an  operation  which  Is  again  highly  dependent  on  the  program 
environment. 


jd 
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3.  Minimal  Storage  Band  Elimination 


The  disadvantages  associated  with  auxiliary  storage  band 
elimination  are  Inherent  In  the  use  of  auxiliary  storage.  In  this 
section,  we  Introduce  a minimal  storage  band  elimination  algorithm  for 
solving  band  systems  of  linear  equations  with  limited  core  storage  and 
NO  auxiliary  storage.  Rather  than  trade  Input/output  and  auxiliary 
storage  for  a reduction  In  core  storage,  we  trade  additional  arithmetic 
2 

operations  Instead.  The  basic  Idea  behind  the  algorithm  Is  quite 
simple.  We  discard  most  nonzero  coefficients  In  the  reduced  triangular 
system  as  they  are  computed  and  regenerate  them  during  back-solution; 
only  enough  Information  Is  retained  to  solve  for  a subset  of  the 
unknowns  at  each  stage,  thus  reducing  the  size  of  the  problem  to  be 
solved. 


Observe  that,  in  solving  (S)  using  auxiliary  storage  band 


elimination,  after  the  (N-m)  variable  has  been  eliminated,  the 
coefficients  and  right  hand  side  of  the  resulting  m x m system  remain  In 
core.  Thus  we  can  solve  for  the  last  m components  of  the  solution 
without  referencing  auxiliary  storage.  Moreover,  we  can  use  these 
values  to  reduce  the  size  of  the  problem  to  be  solved.  Writing  (S)  as  a 
block  system  conforming  to  the  N-m  unknown  variables  x^  and  the  m known 

variables  x^. 


we  see  that  the  remaining  unknowns  x^  satisfy  the  reduced  system 


^11*1 


bi  - Aj2X2  i b^. 


This  suggests  the  following  method  for  solving  band  systems  with 
limited  core  storage  and  no  auxiliary  storage: 

(1)  Use  band  elimination  to  eliminate  the  first  N-m  variables, 

discarding  the  rows  of  U as  they  are  computed  (l.e.,  we  store  only 
the  window) . 


(2)  Solve  the  resulting  dense  positive  definite  system  of  m equations 
in  the  last  ra  variables. 


Advances  in  semiconductor  technology  have  made  the  central  processing 
unit  (as  well  as  attached  processors  such  as  array  processors) 
relatively  Inexpensive  compared  to  auxiliary  storage  devices.  Thus, 
such  an  exchange  Is  not  unreasonable. 


L 


T 


(3)  Substitute  these  values  Into  the  original  systea  to  obtain  a new 
system  of  N-m  equations  In  the  remaining  N-a  unknowns. 

(4)  Solve  the  new  systea  by  the  same  algorithm.^ 


Clearly  the  storage  required  Is  the  same  as  that  for  auxiliary 
storage  band  elimination,  l.e. 

S(N,m)  - 1/2  (n»+-l ) (iiH-2 ) ; 

but  what  Is  the  Increased  cost?  Letting  6(n,m)  denote  the  work 
required,  we  have^ 

e(N,m)  s e (N,m)  + e(N-m,m) 

D 


* 1/2  Nm  + e(N-m,m)  . 


Therefore, 


e(N,m)  a 1/4  N m, 

l.e.,  the  work  has  Increased  by  a factor  of  N/2m.  For  our  model 

3 2 

problem,  we  have  reduced  the  storage  required  from  ~ n to  ' 1/2  n , but 

4 5 

the  work  has  increased  from  ~ 1/2  n to  ” 1/4  n . Fortunately,  we  can 
do  much  better. 

A basic  technique  for  speeding  up  algorithms  is  dlvide-and-conquer . 
Rather  than  reduce  a problem  to  a slightly  smaller  subproblem  at  each 
stage,  it  is  often  more  efficient  to  reduce  it  to  TWO  subproblems  each 
of  at  most  HALF  the  size.  Thus,  suppose  we  could  solve  for  the  middle  m 
components  of  the  solution.  If  we  write  (S)  as  a block  system 
conforming  to  the  first  ~ (N-m)/2  variables  x^,  the  middle  m variables 

x^,  and  the  last  ' (N-m)/2  variables  x^. 


Note  that  the  coefficient  matrix  in  the  new  system,  being  a principal 
submatrix  of  A,  is  again  a symmetric  positive  definite  band  matrix  with 
bandwidth  (at  most)  m;  thus  the  procedure  can  be  applied  to  the  new 
system. 

4 

For  ease  of  exposition,  we  have  assumed  that  the  bandwidth  does  not 
decrease  and  that  only  m components  are  solved  for  at  each  stage; 
otherwise  the  work  would  be  somewhat  smaller. 
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then  the  remaining  unknovms  and  x^  satisfy  two  INDEPENDENT 
'(N-m)/2  X ~(N-m)/2  linear  systems. 


Afi  Xj 

“ **1 

- Aj2  X2 

S 

A33  x^ 

’ '’3 

■ ^32  *2 

3 

^3- 

the  coefficient 

matrices 

of 

the  two  subproblens  are  again 

principal  submatrices  of  A,  they  are  symmetric  positive  definite  band 
matrices  with  bandwidth  (at  most)  m.  Thus  the  algorithm  can  be  applied 
recursively  to  each  of  the  subproblems.  However,  to  use  this  technique, 
we  need  a method  for  solving  for  the  middle  m components  of  the 
solution. 


Recall  that  symmetric  Gaussian  elimination  can  be  viewed  as  an 
elimination  process.  That  Is,  each  equation  is  used  to  eliminate  the 
corresponding  variable  from  the  remaining  equations,  and  the  resulting 
triangular  system  is  solved  for  the  unknown  variables.  For  band 
elimination,  we  observed  that,  when  we  eliminate  the  variables  In 
forward  order,  no  nonzero  coefficients  are  created  outside  the  band. 

Thus  we  need  only  store  and  operate  on  coefficients  within  the  band.  Of 
course,  the  same  Is  true  if  we  eliminate  the  variables  In  reverse  order, 
since  reversing  the  order  of  the  variables  and  equations  preserves  the 
band  structure.  Moreover,  since  A is  positive  definite,  any  order  of 
elimination  Is  numerically  stable  (71. 

This  suggests  another  minimal  storage  band  elimination  method: 

(1)  Eliminate  the  first  "■  (N-m)/2  unknowns  In  forward  order,  using  a 
triangular  window  and  discarding  the  rows  of  U. 

(2)  Eliminate  the  last  "■  (N-m)/2  unknowns  In  reverse  order,  again 
using  a triangular  window.^ 

(3)  Solve  the  resulting  dense  positive  definite  system  of  m equations 
In  the  middle  m unknowns. 


Evans  and  Hatzopoulos  [3] 
solve  centro-symmetrlc  (a^^ 

equations. 


have  used  two-sided  band  elimination  to 

•a  ...  band  systems  of  linear 

n-i+l,n-J+l 


I 
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(4)  Use  these  values  to  split  the  problem  into  tw  Independent 
subproblems. 

(5)  Apply  the  algorithm  recursively  to  solve  each  subproblem. 

The  storage  required  Is  approximately  twice  that  of  our  earlier 
method,  since  we  must  maintain  the  contents  of  both  windows;  in  fact, 

precisely^ 

SMs(N,m)  - (nrt-I)^ 

locations  are  required.  Let  6„-(N,m)  denote  the  %#ork  required  to  solve 

MS 

a system  of  N equations  with  bandwidth  m using  this  algorithm.  The  work 
to  solve  for  the  middle  m unknoms  Is  exactly  the  same  as  before  so 

that^ 

a 0^(N,m)  + 2 ej^g((N-m)/2,m) 

a 1/2  N m^  + 2 . 

Thus 

0Ms(N,m)  a 1/2  Nm^  log2(2N/m). 

While  minimal  storage  band  elimination  avoids  the  problems 
associated  with  the  use  of  auxiliary  storage.  It  does  have  certain 
disadvantages.  First,  the  work  required  has  Increased  logarithmically, 
although  this  can  be  improved  somewhat  by  reordering  the  subproblems  at 
each  stage  to  minimize  bandwidth  (see  Section  4).  Second,  the  method 
has  no  memory:  We  solve  for  one  right  hand  side  but  do  not  save  any  of 
the  Information  needed  to  solve  for  additional  right  hand  sides  (as  In 
auxiliary  storage  band  elimination).  Third,  the  method  is  somewhat  more 
complex  than  In-core  or  atixlllary  storage  band  elimination. 


Storing  two  complete  windows  would  require  (nri-l ) (nrt-2)  locations. 
However,  when  reverse  elimination  (Step  (2))  begins,  the  last  colixnn  of 
the  forward  window  contains  an  unmodified  column  of  A and  need  not  be 
stored. 

^ Again,  for  ease  of  exposition,  we  have  assumed  that  the  bandwidth  does 
not  decrease. 


-J.0- 


4.  Minimal  Storage  Band  Elimination  for  the  Model  Problem 

In  deriving  the  operation  count  for  minimal  storage  band 
elimination,  we  did  not  make  any  additional  assumptions  about  the 
Induced  subproblems.  Ue  merely  observed  that  the  coefficient  matrices, 
being  principal  submatrices  of  A,  were  again  symmetric  and  positive 
definite  with  bandwidth  (at  most)  m.  As  a consequence,  the  decrease  In 
storage  engendered  a logarithmic  Increase  In  work.  However,  If  the 
matrix  A has  more  structure,  one  would  expect  that  the  bandwidths  of 
some  subproblems  could  be  reduced  by  reordering  the  variables  and 
equations,  thus  reducing  the  total  work  required.  The  savings  can  be 
quite  significant.  We  show  In  this  section  that  the  number  of 
arithmetic  operations  needed  to  solve  our  model  problem  Is  approximately 

g 

5/3  that  for  band  elimination. 

Consider  minimal  storage  band  elimination  applied  to  a five-point 
problem  on  an  n x n grid.  During  the  first  partial  solution  step,  we 

9 

solve  for  the  variables  on  a horizontal  grid  line  which  divides  the 
grid  into  two  'n/2  x n subgrids  (see  Figure  4a).  Given  these  values, 
the  proolem  splits  Into  independent  five-point  problems  on  each  of  the 
two  subgrids.  The  optimum  ordering  for  each  subproblem  (in  terms  of 
minimizing  the  bandwidth)  Is  the  column-by-column  ordering  for  which  the 
bandwidth  Is  "nil  [4].  During  the  second  partial  solution  step  (as 
applied  to  the  optimally  ordered  subproblems),  we  solve  for  the 

9 

variables  on  vertical  grid  lines  which  further  subdivide  the  grid  (see 
Figure  4b).  Each  resulting  subproblem  Is  a five-point  problem  on  an 
~n/2  X ~n/2  grid.  Thus  we  have  reduced  our  problem  to  four  structurally 
similar  problems  of  approximately  one-fourth  the  size  which  can  now  be 

1 ^ . w 10 

solved  in  the  same  manner. 


Similar  results  are  easily  obtained  for  more  general  finite  difference 
and  finite  element  approximations  to  self-adjoint  elliptic  partial 
differential  equations.  The  key  property  is  that  the  subprobleros 
generated  have  a sufficiently  regular  structure  that  one  can  analyze  the 
effect  of  reordering  to  minimize  bandwidth. 

9 

When  n Is  even,  the  middle  m variables  do  not  lie  on  a single 
grid-line;  however,  we  could  equally  well  solve  for  the  m variables  on 
either  of  the  two  central  grid-lines. 


It  is  interesting  to  note  that  the  order  in  which  we  solve  for  the  i 
variables  is  essentially  the  reverse  of  the  nested  dissection  ordering  I 
proposed  by  George  [4]  for  sparse  elimination.  j 


-11- 


Let  3 (p)  denote  the  number  of  arithmetic  operations  required  to 

NS 


solve  a five-point  problem  on  a p x p grid  using  optimally  ordered 
minimal  storage  band  elimination.  From  the  preceding  discussion,  we 
have 


eg(p^,p)  + 2 0jj(p^/2,p/2)  + A e^g(p/2) 


1/2  (p^(p)^  + 2(l/2)(p^/2)(p/2)^ 


4 0„s(p/2)  + 5/8  p^  + O(p^). 


Therefore 


5/6  n^  + O(n^) 


so  that  the  total  work  Is  approximately  5/3  that  for  band  elimination. 
The  total  storage  Is  just 


(n+1)^ 


~ 3 

versus  n for  band  elimination. 


During  the  execution  of  this  algorithm,  we  eventually  reach  a value 
of  p sufficiently  small  that  the  entire  factorization  of  each  p x p 
subproblem  can  be  kept  In  core.  In  y.  rtlcular.  If 


p ^ n 


2/3 


then  the  storage  required  for  band  elimination  satisfies 

2 


Sg(p^,p)  = (p^)(p) 


3 ^ , 2/3,3 

p £ (n  ) 


whereas  we  have  assumed  that  the  amount  of  storage  available  Is  at  least 
2 

n . Thus  we  could  switch  to  band  elimination  at  this  point.  The 
switch  would  have  very  little  effect  on  the  total  number  of  arithmetic 
operations  since  most  of  the  work  is  done  in  the  first  few  stages. 
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Flgure  2a.  An  nxn  grid. 


Figure  2b.  The  natural  row-by-row  ordering  of  an  nxn  grid. 

********* 

*1234567* 

* 8 9 10  11  12  13  14  * 

* 15  16  17  18  19  20  21  * 

* 22  23  24  25  26  27  28  * 

* 29  30  31  32  33  34  35  * 

* 36  37  38  39  40  41  42  * 

* 43  44  45  46  47  48  49  * 

********* 
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Pig.  3a.  Symmetric  band  elimination. 


1 


I 


FOR  k - I UNTIL  N DO 
{ Imax  ■ MIN  (k+m,  N); 

FOR  1 - k UNTIL  Imax  DO 
M[k,ll  - Alk.l] ; 
x(kl  - b[k] }; 


FOR  k - 1 UNTIL  N DO 
{Imax  “ MIN  (k+m,  N); 

FOR  1 - k+-l  UNTIL  Imax  DO 
{ulk  - M[k,l]  / M[k,kl ; 

FOR  j - 1 UNTIL  Imax  DO 

Md.JI  - Md.Jl  - ulk  * Mlk.jl  ; 
x(ll  - x(l]  - ulk  * x[k) ; 

M(k,ll  - ulk); 
x[k]  - x[k]  / M[k,k] ); 

FOR  k - N STEP  -1  UNTIL  1 DO 
{imax  - MIN  (k+m,  N); 

FOR  1 - k+-l  UNTIL  Imax  DO 

xtk]  - x(k)  - M(k,l]  * xdl  }; 


Figure  3b.  The  window  on  the  band  of  A advances  at  each  step. 
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Figure  4a.  The  variables  on  a horizontal  grid  line  subdivide  the  grid. 
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Figure  4b.  The  variables  on  vertical  grid  lines  further  subdivide  the  grid.  j 
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